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An analog computer was used to simulate zone melting in small-diameter tungsten rods 
of 0.2 to 1 centimeter. Heat balance equations take account of electron bombardment, heat 
capacity, radial conduction, axial conduction, and surface radiation. Diagrams are given 
for the chosen spatial subdivision and for the computer mechanization. Predicted 
temperature-time curves are shown for a rod of 0.4-centimeter diameter with optimum 
power input of 650 watts. Unsatisfactory operation is predicted for zone melting of a 
1-centimeter rod. 



As part of the NBS program for attempting- to 
supply metal samples so pure that no spurious lines 
can be detected by spectroscopy, G. A. Moore of 
the Chemical Metallurgy Section has developed 
an apparatus for zone melting of small-diameter 
tungsten rods. Because of the high cost of the semi- 
purilied material used in the process, it became 
desirable to simulate the heat transfer operation as 
a means of determining optimum relations between 
rod diameter, heating rate, and progress of the molten 
zone. 

The technique is applicable to other materials, 
and the problem is interesting from the simulation 
viewpoint for several reasons. The conventional 
diffusion equation frequently used for analog com- 
putation is not suitable because it requires precise 
knowledge of the temperature distribution at the 
surface. A model must be set up in terms of heat 
flows; these too are unknown, but can more readily 
be estimated from the physical situation. The 
model must include radiation losses, which at these 
high temperatures require fourth-power nonlinear 
terms. The change of state at the melting point 
results in a local interruption of the temperature 
rise in the melting zone. Conductivity," thermal 
capacity, and emissivity will change with temperature 
and with state. Time rates of temperature change 
will be treated as infinitesimals, spatial rates as 
discrete increments. The incremental dimensions 
for radial and axial subdivision must be carefully 
chosen if optimum use is to be made of the limited 
number of integrating amplifiers available. 

The radii r s of the 5 rods to be simulated range from 
0.1 to 0.5 cm; lengths (actually 3 or 4 ft) are assumed 
infinite. Heating is by continuous electron bom- 
bardment in a vacuum from a heated wire of 1 .2 mm 
diameter, shaped to an annulus of 4.0-cm diameter 
surrounding the rod as shown in figure 1. The 
heater is initially fixed ' near one end, but after 
melting begins it is moved very slowly up (or down) 
the rod. The impurities are retained in the molten 
zone and are carried to one end, leaving the purer 
metal to solidify again to rod form. r I his motion 
is neglected in the simulation. Circumferential 
symmetry and an isomorphic rod are assumed, and 
the heat input from the wire is taken to be a step 
function of time. 



Because of symmetry on either side of the heater 
midplane, it is enough to consider only half the length 
of the rod. Preliminary estimates indicated that 
temperatures would rise only slightly at points 20 mm 
distant from that plane, so that heat flow beyond this 
was neglected. Steep temperature gradients may be 
expected near the heater, and it is necessary to assure 
that metal will not drip off the outer surface before 
the core is melted. One purpose of the simulation, 
therefore, was to find whether the molten zone ad- 
vanced radially at least as fast as it advanced axially. 
Subdivisions into annuli or cylinders with lengths A?/, 
outer radii r , and inner radii r t were decided upon. 
Their generalized proportions relative to surface radii 
r s are given in table 1 . 




>> HEATER WIRE 



Figure 1. Assumed subdivision of rod. 
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The heat stored per unit time in zone 1, which is 
assumed to have a uniform temperature of T\ degrees 
Kelvin, at time 6 seconds after heating begins, is the 
thermal capacity of the zone volume multiplied by 
the rate of temperature rise. With density yg/cm 3 
and temperature-dependent specific heat c=c(Ti) 
cal/g-deg K, this storage rate is ycir(;r 2 s —rV)Ay dTijdB. 
Similar expressions hold for the other zones, with r 
replacing r s and r*=0 for the inner ones. 





Table 1. 




Zone 


r /r s 


Ti/Ts 


Ay 








Mm 


1 


1 


0.5 


1 


2 


0.5 





1 


3 


1 


0.5 


1 


4 


0.5 





1 


5 


1 


0.5 


3 


6 


0.5 





3 


7 


1 





5 


8 


1 





5 


9 


1 
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Heat inflow to zone 1 will be 239 P cal/sec, where P 
is the steady kilowatt power input to the half-rod. 
Heat outflow by reason of conductivity k=k(T 1 ) cal/ 
(sec) (cm) 2 (deg K/cm) will be radially inward to 
zone 2 across the area 2Trr f Ay with the gradient (T Y — 
T 2 )/(ri — r 2 ), where n and r 2 are mean radii. Also, 
there will be conduction axially to zone 3, at the rate 
wk(r 2 8 —r() {T x — T z )l{y x —y z ) where y x and y 3 are mean 
lengths. No heat is transferred beyond the heater 
midplane, since there is no gradient across that 
boundary. However, zone 1 loses heat by radiation 
to ambient temperature T a at the rate 1.37e (2irr s Ay) 
[(T 1 /1000) 4 -(T G /1000) 4 ] where e=e(T 1 ) is the tem- 
perature-dependent emissivity. The heat flow bal- 
ance equation for zone 1, at any time up to the 
instant of local melting, is then 



7 ar(r, 2 -r?)Ay§=239P 



- 27Tr 7 i 



--*«-'•> [Ml 

-1.37.(2^!,) [(jgu)'-^)'} 

Similar heat flow balance equations can be written 
for the other zones except that 239P is replaced by 
radial or axial conductive inputs, and there is no 
radiation for the interior zones where r replaces r s and 

r«==0. 

For zone 3, for example, the equation is 

ycr(r*-rVAy^=-2^y[^] 

-1.37e(2 ir r s A 2/ )[( I ^) 4 -( I §- ) 4 ] 



Numerical values of the properties y, c, k, and e 
were known only imprecisely for tungsten, but it ap- 
peared that they change only slowly in the range 
from T a =300 to 7^=3700, the melting point. Ac- 
cordingly, the values of table 2 were used in the com- 
putation, based on preliminary estimates of temper- 
atures with zone 1 melted and zones 2 and 3 just 
beginning to melt. The heat of fusion in zone 1, 
about 44 cal/g, was neglected in this approximation. 
For each zone, the initial temperature was assumed 
to be 300 °K, and regions beyond zone 9 were as- 
sumed to remain at this temperature. 

Substitution of numerical values in the 9 expres- 
sions for heat flow gave for the rod of 1-cm diameter 
the following set of simultaneous first-order differ- 
ence-differential equations for solution by the analog 
computer : 

+ 0.2T 1 = .015P-.0394(r L -7 7 2 )-.369(T 1 -T 3 ) 
-.0498T?/10 6 

-2T 2 =-.118(r i -T 2 ) + .363(r 2 -r 4 ) 
-0.2T 3 = -.369(7!- 7 7 3 ) + .0388(T3-T 4 ) + .173(r 3 -T 5 ) 

+ .0498TI/10 6 
+ 0.2T 4 = + .406(T 2 -r 4 ) + .1302(T 3 -r 4 )-.23(J r 4 -r 6 ) 
+ 0.2T 5 - + .439(2 T 3 -T 5 )-.3195(T 5 -T 6 )-.261(r 5 -r 7 ) 
-.134rt/10 6 
-2T 6 =-.115(T 4 -T 6 )-.2116(T 5 -2 7 6 ) + .067(2 7 6 -r 7 ) 
-T 7 =-.129(r 5 -T 7 )-.0504(T 6 -T 7 ) + .167(2 7 7 -r 8 ) 
-f= -.W7(T 7 -T 8 ) + AQ7(T S -T 9 ) 
+ T 9 = + .167(T 8 -T 9 )-.167T 9 + .502 



Table 2 



Zone 


T 


7 c 


k 


e 


1 


3700 


0.84 


0.31 


0.36 


2 


3700 


.84 


.31 




3 


3700 


.84 


.31 


.36 


4 


3000 


.75 


.30 




5 


1100 


.55 


.27 


.13 


6 


700 


.50 


.39 




7,8,9 


500 


.50 


.42 


.05 



Figure 2 shows the computer patchboard diagram 
for the first two and last two of these equations. 
The time scale was slowed in the ratio of 1:20, and 
1 machine volt represented 100° Kelvin. Heater 
input P was adjusted by trial with the aim of estab- 
lishing an optimum value which would use a mini- 
mum of energy without causing zone 3 to melt 
ahead of zone 2. For the 1-cm rod, such a value 
could not be found, since zone 3 always reached 
3,700° before zone 2, while the temperature of zone 
1 went up to about 4,600°. This represents an 
undesirable physical situation which in practice may 
prevent satisfactory processing of a rod of this large 
diameter. In view of the uncertainty as to actual 
values of y, c, k, and e, it was considered unnecessary 
to carry through the computation for a more refined 
approximation. 
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Figure 2. Initial and final stages of computer patching. 



A sot of temperature-time curves for one of the 
smaller rods with a 650-watt input is shown in 
figure 3. It will be noted that the melting of /one 
3 lag's that of zones 1 and 2, and hence that the 
power can be adjusted so that zone 3 just fails to 
melt. This assures formation of a thin liquid layer 
which is physically stable. The rod is nearly melted 
through in about 3 sec. This finding required modi- 
fication of the power control circuits of the experi- 




5 6 

TIME ,sec 



Figure 3. Zone melting 2-mm-radius tungsten rod. 

mental apparatus by addition of rapid-response 
components to supplement the original servo- 
mechanism, which had required about 30 sec to 
establish an equilibrium setting. 



Assistance on this problem came from NBS 
colleagues G. A. Moore, who supplied the physical 
constants and offered much constructive advice; P. 
K. Wong, who operated the computer; and K. 
Shives, who plotted the curves. 

(Paper 65C2-60) 
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